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Abstract. First principles simulations of the quantum dynamics of interacting Bose 
gases using the stochastic gauge representation are analyzed. In a companion paper, 
we showed how the positive P representation can be applied to these problems using 
stochastic differential equations. That method, however, is limited by increased 
sampling error as time evolves. Here, we show how the sampling error can be greatly 
reduced and the simulation time significantly extended using stochastic gauges. In 
particular, local stochastic gauges (a subset) are investigated. Improvements are 
confirmed in numerical calculations of single-, double-, and multi-mode systems in 
the weak mode coupling regime. Convergence issues are investigated, including 
the recognition of two modes by which stochastic equations produced by phase 
space methods in general can diverge: movable singularities and a noise-weight 
relationship. The example calculated here displays wave-like behaviour in spatial 
correlation functions propagating in a uniform ID gas after a sudden change in the 
coupling constant. This could in principle be tested experimentally using Feshbach 
resonance methods. 
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1. Introduction 

Phase-space techniques El for the simulation of exact many-body quantum dynamics 
have been developing rapidly^ El El El EI- Perhaps the major driving force behind this 
is their ability to work around the exponential growth§ of Hilbert space that occurs in 
direct methods that work with an explicit quantum stated • This exponential scaling in 
the direct approach implies that it is virtually impossible to diagonalize general initial 
states in order to calculate dynamics, even in the exceptional case when the exact 
eigenstates are known. Path integral methods were developed to deal with the Hilbert 
space growth, but while they are very useful for calculating ground states (TUIIIII, they 

I www.physics.uq.edu.au/BEC 

§ Exponential in the number of modes or orbitals involved. 
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are not applicable to quantum dynamics due to the oscillatory phase terms which arise 
and swamp any observable predictions. In contrast, appropriate phase-space methods 
can scale linearly or polynomially with the system size (like path integrals) without being 
as affected by the oscillatory phases, and so can provide verifiable physical predictions. 

In the companion paper ^2], we analyzed the performance of the nonclassical 
positive-P representation in U] for many-body interacting bosonic systems. Here, we 
extend the analysis to treat the more flexible stochastic gauge method of weighted 
phase-space trajectoriesjZj. 

The aforementioned positive P representation is a straightforward and successful 
method, which is a good "baseline" on which more involved methods can build. It is 
useful to give a brief overview of it and its limitations, at this point: The quantum state 
is written as a probability distribution of off-diagonal coherent states separable at each 
spatial lattice point. This corresponds to a nonclassical phase-space with twice the usual 
classical dimension. The number of variables needed to specify a single coherent state 
configuration grows only linearly with the lattice size because they are separable, while 
correlations between subsystems are contained in the details of the distribution (which 
is stochastically sampled). This very mild growth of the number of required variables is 
the reason that this method can lead to tractable many-body simulations. 

Dynamics in the positive-P representation takes on the particularly simple form 
of the widely used Gross-Pitaevskii mean-field equations [in], but with appropriate 
Gaussian noise terms added. Several interacting Bose gas dynamics simulations with 
this method have made quantitative predictions, including: 1) The scattering dynamics 
during the collision of two BECs, and the evolution of momentum correlations between 
the scattered atoms [T6]. 2) Dynamics of evaporative cooling and incipient condensation 
of an interacting Bose gas^lE]; 3) Dynamical behaviour of spatial correlations in one- 
dimensional interacting Bose gases HE] : 4) Dynamics of quantum soliton propagation 
in optical fibers pT ] IT ^ I19j : and 5) Dynamics and steady-state behaviour of open boson 
systems coupled to reservoirs [B]; Dynamics and thermodynamic spatial correlations of 
a similar but explicitly particle-conserving model have also been simulated using the 
related stochastic wavefunction method [5] 1^. 

While much can be simulated with the positive P, there is a limiting factor — the 
growth of the scatter of trajectories with time. This statistical error can be estimated, 
but eventually after some evolution time it becomes large enough to obscure the physical 
behaviour for any practical number of trajectories. Estimates of the useful simulation 
time with the positive P method have been obtained^^], and indicate two major 
limitations: 

(i) When occupation of individual lattice points is high 1), the simulation loses 
precision before phase-collapse can occur. 

(ii) Useful simulation times decrease with reduced lattice spacing. 

This loss of precision is often (but not always) associated with the development of power- 
law tails in the distribution function, which can also lead to systematic boundary term 
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errors. A typical cause of power-law tails is the presence of movable singularities in 
the dimension-doubled phase-space equations (discussed in Section IH)). These problems 
are intimately related, and therefore the methods that improve precision can also solve 
boundary term issues. 

One approach to this problem is to change the basis. For example, a different 
phase-space technique that can have lower sampling error when mode occupations are 
high is the truncated Wigner representation^ QHl sometimes called the classical 
field method|22j. The simulation can be quite accurate, particularly when all modes are 
strongly occupied. However, this method involves truncation of third-order differential 
terms in the evolution. This can lead to large systematic errors which are essentially 
intrinsic to the method itself, and tend to grow with time The conditions under 
which it is valid for interacting Bose gases have been investigated in significant detail 
in Ref. 

Here, however, we will investigate a different avenue of improvement. For 
a given basis, there is a wide range of freedom in the correspondence between 
quantum mechanics and the stochastic equations[71 |H]. The available freedoms can 
be constructively specified in a unified way by the stochastic gauge formalism[3 ITHj . 
which adds freely defined gauge functions and a corresponding global weight to the set 
of stochastic variables. Preliminary attempts at harnessing this freedom in single-mode 
systems have been successful j2 5 1 Ej, with demonstrations of time-reversible quantum 
simulations with up to 10^^ interacting particles J26j. In this paper we investigate 
these methods systematically, and also present the progress that has been achieved 
in exploiting stochastic gauges for large multi-mode systems. 

We restrict ourselves to the case of gauge functions defined locally (i.e. separately 
for each lattice point), which is the simplest useful case. Sections El and El investigate 
two gauge approaches, each with its own merits. Subsequently, in Sections El and |H1 
the performance of the proposed improved methods is assessed for single-mode, double- 
mode, and multi-mode cases, respectively. Section (HI investigates convergence issues and 
identifies two ways by which stochastic equations produced by phase space methods in 
general can diverge: movable singularities and a noise-weight relationship. The last 
Section simulates the uniform ID gas as a non-trivial example. In particular, the 
dynamics of second order spatial correlations g^'^^x) after a sudden increase in the 
coupling constant at t = are calculated as in the companion paper [T^. 

2. Model 

2.1. The lattice model 

We wish to solve for the quantum time-evolution of a dilute interacting Bose gas. 
Following the companion paper J2|; we consider a lattice Hamiltonian that contains all 
the essential features of a continuum model, provided the lattice spacing is sufficiently 
small. For a rarefied gas of the kind occurring in contemporary BEG experiments, s- 
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wave scattering dominates j27], and the s-wave scattering length is much smaller than 
all other relevant length scales. If the lattice spacing is also much larger than a^, then 
the two-body scattering is well described by an interaction local at each lattice point. 
Otherwise, a more careful renormalization procedure |28j than the one below is required. 

Let us label the spatial dimensions by c? = 1, . . . , P, and label lattice points by the 
vectors n = (ni, . . . , ud). For lattice spacings Ax^, the spatial coordinates of the lattice 
points are Xn = (niAxi, . . . ,n£iAx£)). The volume per lattice point is AV = Ax^. 
This lattice implies effective momentum cutoffs[2ni of fc™'^^ = tt/Ax^. We also define 
the lattice annihilation operators (~ V AV ^E'(xn) in the field notation of Eqn. (1) of 
[12]), which obey the usual boson commutation relations of [an,aL] = "^nm- With these 
definitions, one obtains: 



H = h 



(1) 



In this Hamiltonian, the frequency terms c^nm = ^mn come from the kinetic energy 
and external potential. They produce a local particle-number dependent energy, and 
linear coupling to other sites, the latter arising only from the kinetic processes. The 
nonlinearity due to local particle-particle interactions is of strength 

with the standard coupling value p7] being g = g^j^ = Anash^ /m in 3D, and g = g^D/cr 
in 2D and ID, where a is the effective thickness or cross-section of the collapsed 
dimensions. 

When interaction with the environment is Markovian (i.e. no feedback), the 
evolution of the density matrix p can be written as a master equation jHOl 1^ in 
Lindblad^ form 

dp 1 



dt ih 



For example, single-particle losses at rate 7n (at Xn) to a T = heat bath are described 
by Ln = a^y/j^. 

2.2. Gauge P representation 

The gauge P representation was introduced in Ref. jH] and is described in more detail 
in Refs. [ZllEl- Here we summarize the issues relevant to the dynamics of the model 
(H)), and present the stochastic equations to simulate. 

For a lattice with M points, the density matrix is expanded as 

p = j G{a,P,n)A{a,/3,n)£^^a£^^f3£n, (4) 

where A is an off-diagonal operator kernel, separable between the M lattice point 
subsystems. We use a coherent state basis so that 

A{a,f3,n) =^]||a)(/3*||exp[-a■/3], (5) 
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in terms of Bargmann coherent states with complex amphtudes o: = ( 

1 1 «) = ®n exp [anol] I 0) . (6) 

n is a complex weight. 

The gauge representation is a generalization of the positive P representation, which 
has no weight term. With the choice 



(7) 



we recover the positive P distribution P+. It has been shown that all density matrices 
can be represented by a positive real -P+jS], so the same is true for G. The expansion 
(13)) then becomes a probability distribution of the A, or equivalently the variables a, 
f3, fl. A constructive expression for a distribution P+{p) is given by expression (3.7) in 
Ref. 0, although more compact distributions may exist. In particular, a coherent state 
\a.o) will simply have 



G = 6''' (a - cto) 6'"' {(3 - a*,) 6\n - 1). 



Using the identities 







/3n + 



n 



_d_ 

d 



A, 



A, 



(8) 

(9a) 
{9b) 

(9c) 



the master equation © in p can be shown to be equivalent to a Fokker Planck equation 
in G, and then to stochastic equations in the an, (i^i ^^^1 f2, provided that boundary 
terms vanish in carrying out the partial integration used to derive the Fokker-Planck 
equation. This implies a restriction on the phase-space distribution tails, which typically 
must vanish faster than any power law. 

The standard method is described in Ref. ^Slj. The correspondence is in the 
sense that appropriate stochastic averages of these variables correspond to quantum 
expectation values in the limit when the number of trajectories S oo. In particular 
one finds thatjO] 



lim 



(10) 



Any observable can be written as a linear combination of terms (jlOj) . We will use the 
notation (■)s to distinguish averages of random variables from quantum expectations 

(■)• 

Due to the fact that the basis vectors \\a.) are not mutually orthogonal, many 
different distributions (and hence, sets of equations) correspond to the same master 
equation. It has been foundfUIinj that the family of stochastic equations corresponding 
to a given master equation is parameterized by stochastic gauge functions of several 
kinds. These are completely arbitrary functions that appear in the equations, but do 
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not affect However, the rate at which precision of the estimates (|iup improves with 
more stochastic reahzations can be strongly affected. It is shown below that judicious 
choices of the gauges can lead to large improvements in precision. 

The 2M + 1 Ito stochastic equations to simulate are found using methods described 
in Refs. jHH Ej . For the model ((H) obeying the master equation © with coupling to a 
zero temperature heat bath, they are 

dan = [-i ^ t^nmam - inUnan - 7nan/2](it + ^ B^^^{dWk - Qkdt), (11a) 



k 



dPn = Xl<m/?m + ^«n„/?„ - 7n/5„/2]rft + ^ilV^^fc " Okdt), (116) 



dQ= n 



Y^GkdWk 



[lie) 

k 

Here, = anPn, and there are M' > 2M labels k for the independent noise 
terms, to sum over. The dWk are independent Wiener increments {dWj{t) dWk{s))s = 
5jkS{t — s)dt^ . In practice, these can be realized at each time step At by mutually- and 
time-independent real Gaussian noises of mean zero, and of variance At. The elements 
of the M X M' noise matrices B must satisfy 

k 

E^S^^mi =^'^/3Xm, (126) 
k 

E^iM=0. (12c) 

k 

The M complex quantities Qk are arbitrary complex gauge functions, referred to 
as drift stochastic gauges. There is also more gauge freedom here because ()2.2j) do not 
specify the noise matrices i?^"-* and B^^^ precisely. This freedom can be expressed as 
diffusion stochastic ga,uges^I^W\Yi\. 

The simulation strategy is (briefly): 

(i) Sample a trajectory according to the known initial condition G(0) = (j'(p(0) ) 

(ii) Evolve according to the stochastic equations ()2.2|) . calculating moments of interest, 
and recording. 

(iii) Repeat for 5 ^ 1 independent trajectories and average. 
2.3. Single-mode model 

The single mode model is a good test bed for gauge choices in quantum many- 
body theory. It is exactly solvable, yet it already contains all the essential features 
that lead to the rapid growth of fluctuations which limits simulation time with the 
positive P method^2]- It has also been experimentally realized in recent optical lattice 
experiments |33j with Bose-Einstein condensates. We proceed as in the companion 
paper [T2j: 
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To simplify the notation, the mode labels n = (0, . . . , 0) will be omitted when 
referring to the single mode system. Furthermore, we move to an interaction picture 
where the harmonic oscillator evolution due to the ua^a term in the Hamiltonian 
is implicitly contained within the Heisenberg evolution of the operators. Then, this 
"anharmonic oscillator" model simply has 

H = ^a^^a\ (13) 

When dealing with this system it will be convenient to consider the evolution 
starting with an off-diagonal coherent-state kernel 

Ao = A(ao,/5o,l), (14) 

with "particle number" no = ao/3o and initial unit weight i7(0) = 1. This is because 
for any general initial state, each sampled trajectory will start out as an Aq with some 
coherent amplitudes. 

With initial condition Aq, analytic expressions for observables can be readily 
obtained. In particular, the first-order time-correlation function 

G«(0,t)=/5o(a)=a*(at)*, (15) 
which contains phase coherence information. Normalizing by Tr[Ao], 

G«(0,t) = noe-^*/2 / (g— <-7* _ i) \ (ig) 

1 1 - n/t J 
When the damping is negligible, no real, and the number of particles is no ^ 1, 
one sees that the initial phase oscillation period is 

tosc = -sm"M — , (17) 

and the phase collapse timejJJ over which [(^^^^(O,^)! decaysf is 

When there is no damping, the first quantum revival occurs - as observed 
experimentally inn] - at 

revival ~ ~' 

(19) 

2.4- Positive P noise behaviour 

The noise behaviour of this model and simulation times has been investigated in detail 
for the positive P case in Ref. [T2j. Some useful conclusions include: 

f We have kept the notation icoh from the companion paper [T^. where "coh" indicates that coherence 
is maintained. 
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(i) If ^ is a Gaussian random variable, and v = e"^, then finite-sample estimates of (f )g 
using S realizations of the random variable have relative fluctuations that scale as 



(ii) In a single mode system, the variable n{t) is the exponential of a Gaussian random 
variable, as is a{t) at short times. Since observable estimates (fTUI) involve means 
of polynomial functions / of these variables (/)s; then one needs 



for reasonable precision. 

(iii) The mechanism responsible for limiting simulation time is growth of fluctuations 
in loglal due to the real part of the d\oga = —inndt term. When \n\ ^ 1, the 
fluctuations in Im [n] are sizeable, and var [log |a|] rapidly reaches large values that 
exceed (1^ . 

(iv) In the multi-mode system, this nonlinear (in a^) term becomes da^ = —inrinan dt, 
which depends only on the variables in the local mode n. For this reason, the 
single mode analysis continues to give a good qualitative description of the noise 
behaviour in multi-mode systems, especially while the inter-mode coupling is weak 
in comparison with the local scattering. 



In what follows, an estimate tsim of the useful simulation time will be evaluated using the 
condition (|7H) for the amplitude variables a and /?. For the positive P method applied 
to the system (fT3|) . tgim was found to be relatively shortest at high mode occupation 
and weak damping, scaling as oc n~^/^. This is shown in Figure 01 and Table H For 
^ 1 we see that the positive P method does not reach even the phase collapse times 
^coh|2II- This is caused by the phase- diffusion inherent in any initial state with a range 
of particle numbers, where the self-interactions cause the phase to evolve at a different 
rate for different total particle number. 

However, indications that simulation improvements can be obtained using nonzero 
gauge choices were seen by Deuar and Drummondj34j and Plimak et al [23], using drift 
and diffusion modifications, respectively. Carusotto et a/[Hj also found improvements in 
a different (number-conserving) model, using an approach that effectively uses a drift 
gauge (see Section E31 for more on this). 

The aim here is to develop gauges that improve simulation times at high mode 
occupation, while being applicable to the many-mode situation where conditions for a 
single mode are dynamically changing. To do this we wish to find the dependence of 
advantageous gauges on mode parameters such as n{t), constants k, 7, and time t itself. 
We are especially interested in the case of low or absent damping 7 <^ k, in highly 





(21) 



2.5. Aims 
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occupied modes n ^ 1. This is the regime where the simulation time is most hmited 
with the positive P method. 

3. Local diffusion gauges 

In the positive P simulation of modes with n ^ 1 occupations, tsim is limited by the 
following process: In order to allow for an exact simulation of quantum fluctuations, the 
nonclassical phase space requires the effective particle number n = n' + in" to develop 
complex values. This in turn gives rise to exponential growth in amplitude fluctuations 
in either |a| or since the drift equations have the structure: 
d\a\ 



dt 

d\P\ 



Kn"\a\ + . . . (22a) 
-Kn"\(3\ + ... (226) 



dt 

These instabilities lead to growth in sampling error, and possibly ultimately to 
systematic boundary term errors as well. In this section we consider the freedoms 
present in the choice of noise matrix B, to limit the stochastic growth in the effective 
'gain' n", without introducing any drift gauges. Such techniques can extend the useful 
time-scale of a simulation, although they are typically unable to remove boundary terms 
at long times if caused by movable singularities. 

In this case = 0, dfl = 0, so there is no weight evolution and we can assume 
f2 = 1 for all practical purposes. Analytic expressions for an optimized diffusion gauge 
choice will be found and discussed. Their performance in actual numerical simulations 
is reported in Sections El to |H1 

3.1. Diffusion gauge mechanism 

The noise matrices appearing in the stochastic equations must obey ()2.2|) which can be 
written matrix equation BB^ = D, where 



B 



(23) 



and D is completely determined by the master equation. Otherwise the B are free, 
including the number of noises M', i.e. columns in the matrix. These freedoms in 
B choices are investigated in Ref. 116J. A broad class of suitable B which are square 
complex orthogonal matrices have been described in These are given by 

B{g,u) = ^0{g,k), (24) 

where O are arbitrary complex 2M x 2M orthogonal matrices such that 00^ = I. In 
the general case, these can be constructed using (2M — 1)M complex diffusion gauge 
functions gjk- 



= exp , (25) 
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a 



Ip 



10 



(26) 



3.2. Single mode 

We proceed similarly to the companion paper^] for the positive P case. For the model 
(fT^ with damping to a zero temperature heat bath, the equations to simulate are 

da = —a[iKn + 7/2] dt + iVina [cos{g)dWi + sm{g)dW2\ 

d(3 = (3[iKn - 7/2] dt + Vin f3 [- sm{g)dWi + cos{g)dW2]. 

with n = af3. Here we include only one complex diffusion gauge g = g^ = g' + ig" ■ 



(27) 



3.2.1. Real gauges Note that 

B = BoO{g)=BoS{g")R{g'). 



(2^ 



in terms of a rotation R and a transformation 5*. The rotation serves only to mix the 
noises dWj together: 



R{g')dW 



cos (7' 


sin g' 




dWi 




dWi 


— sin g' 


cos (7' 




dW2 




dW2 



(29) 



where the new noises dWj have the same statistical properties as the original noises 
dWj. Because of this, g' has no impact on the statistical properties of the simulation. 
Henceforth we will consider only g = ig", where cosg = cosh^f", and sin^f = isinhg". 



3.2.2. Logarithmic variance As in ^2] for the standard positive P case, we consider 
the mean phase variable variance 
1 



V = - {var [log \a{t)\] + var [log |/3(t)|]} 



(30) 



which was found to be the limiting factor for simulation time when it exceeds V > O (10) 
as per (ET| . 

Taking g = ig" the equations (jTfjl can be formally solved using the rules of Ito 
calculus as 



logn(t) =logn(0) -7t + Vme-^ [C2(t) + <i(t)] , (31a) 

loga(t) = loga(O) + (m — ■y)t/2 + iViK(^i(t) cosh (7" — \fiKC,2{t) sinh^f" — in j n{s) ds, 

Jo 

(316) 



where 

Q{t) = f dW{s) 

Js=Q 

are Gaussian-distributed random variables of mean zero, and: 
(Ci(^)Cfc(-s))s = Sjkmm [t, s] . 



(31c) 
(32) 
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It is convenient to define 

no = ao/3o = n'g + iriQ. (33) 

Note that n{t) is the exponential of a Gaussian random variable, and that if a generic 
variable ^ is Gaussian then 

lim (m"')s = — 7T7-^ (34) 

for even m, and zero for odd m. From this, \ims-,oo{G^^)s = exp(cr^/2) and 
lim5_,oo(cos(cr^))s = exp(— cr^/2). Using these we obtain that the limiting variance 

V in the limit of a large number of samples 5 ^ oo is given by: 

V = Ktcosh{2g")/2 - /tVp I ^ ~ ^^i^ ^ ^'^^ I (35a) 

?7 I (? - 7) 27 ^ 

where 

g = 2(7 - /te"^^") (356) 

is a "damping strength" parameter. 

At short times the first term (direct fluctuations from noise in (ilog|a| etc.) 
certainly dominates. For long time-scales the indirect noise mediated by a growing 
spread in n" (remaining terms) has the most serious consequences, since it causes an 
exponential growth in the sampled fluctuations. 

3.2.3. Optimum gauge The first term in ()35aj) (due to direct noise in log a and log/3) 
grows with while the later terms decrease. This indicates that there is a trade-off 
parameterized by g" . The optimum choice is given by a balance between the direct noise 
in the amplitudes, and the indirect noise in the gain term n" . The lowest fluctuations 

V at a given time t will occur when dV{t)/dg" = 0. To choose g", we must decide upon 
a "target time" t = ^Qp^ at which to minimize V. 

When we aim for relatively short target times t^-p^ such that l^l^opt ^ ^ 
2KtQp^e~^^" <^ 1, the optimized gauge is approximated by 

(2Ktopt|nop2 



9opt ~ 5'approx - ^ log 

where, on defining v = ■ytopt, the coefficient a2{v) is 

3 /e^2^-(^3^2t;) + l-4e'^ 



■«2(7iopt) + 1 



(36a) 



a,iv) = - ^-^ ^ j . (366) 

This reduces to 02(0) = 1 in the undamped case. The discrepancy between (36) and the 
exact optimization g'^-p^- found using dV/dg" = is shown for real uq = in FigureH] It 
can be seen that for Uq > O (10) and/or for times shorter than singly-occupied coherence 
time the approximate expression for the optimized gauge choice is still useful. 
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Figure 1. Discrepancy A = g'^^^ — 5approx between g'^^^ (the exact optimization 
of g" for the diffusion-gauge-only case) and the approximate expression (36) for an 
undamped mode with diagonal coherent initial occupation no = tiq. Discrepancy 
values A are shown as contours. 



3.2.4- Useful simulation times By inspection of (31), the behaviour of log |a| and log \j3\ 
is Gaussian-like due to the Q terms. Since G^^^ oc {a)s, the condition (pT|) implies that 
useful precision in G^^^ is obtainable only while V < 10. 

In the simplest case of no damping and coherent state initial conditions at large 
occupation hq = u'q 1, for target times times topt ^ 3tosc/4vr one has 6"^^"??"^°" ~ 
\/3/2Ktoptno- These are times longer than an oscillation period but shorter than 
coherence time tcoh- Taking then the terms of V{t,g") from (36) at t = topt of highest 
order in no, and imposing the precision limit V(tsim5 fi'approx) ~ one obtains the 
expected simulation time for an undamped system as 

W^O(10)W (37) 
This is a large improvement at high mode occupations — compare to the positive P 
{g" = 0) results in Tabled 

Checking back, we see that when the target time is given by topt = ^sim from (jSTj), 
it is indeed much longer than an oscillation period, and also that the terms of V highest 

— 1/2 

order in uq remain so after taking /tt < oc into account. The clearest evidence for 

the validity of the above reasoning — which involves a few approximations — is that 
the numerical simulations of Section d agree with ()37|) . 

At low occupations and with uq = Uq, on the other hand, (^approx ~^ 0; ^^^1 one 
again expects the same simulation time as with the standard positive P equations. 

3.2.5. Strong damping If damping is present, then with a large enough gauge 
g" > |log(K/7) such that g > the regime of linear increase in sampling variance 
can always be reached for times gt ^ 1. Here, fluctuations grow slowly as 

V = Ktcosh(2/)/2-6 . (38) 
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The constant b is 



2n'o - |no|V(e'^" 



(39) 



where e = K/7 < e^^ ". The required g > imphes 7 > Ke~'^^" , i.e. either that damping 
rates are strong compared to the nonhnear detuning at the two-particle level, or the 
diffusion gauge is large. 



3.3. Extension to many-modes 

The expression (36) was worked out under the assumption that the mean occupation of 
the mode is conserved. In coupled-mode simulations this is no longer the case, and can 
be adapted for by replacing hq in the expression (36) by nn(t). 

Next, consider the situation when we aim to simulate until a target time topt- 
some intermediate time < t < t^^^ it may be advantageous to optimize g" only for 
the remaining time to target 

trem = max[topt - i, 0], (40) 

rather than for a constant target time ^Qp^ ahead of t. Indeed, this is demonstrated to 
be an improvement in Section 15.31 

There will be a separate local gauge at each lattice point, such that for an M-mode 
system (with j = 1, . . . , M labeling the modes) the only nonzero noise matrix elements 
are 

iy/inaaj cosh^f^^^ 
-V^anj sinh5f^_^. 
iy/Ji^Pn^ sinh (yf^_^, 
ViKpn, coshg'^^. 
The suggested diffusion gauges are 

■(2/€trem|ri„(t)|)' 



^n,,(i+M) 

r(/3) 
n,,{i+M) 



(41) 



9n 



rem) 



(42) 



3.4. Effect of lattice spacing 

Labeling a or /3 as z, it was found [T^ that for a uniform gas of density p with volume 
AV per lattice point in a D-dimensional lattice one has 

>kinetic^^^^)|^ 



var 



var 



|5direct^^(i)| 



7^ ^ o 



(43) 



at short times (such that the first term in ()35a|l dominates) when using the standard 
positive P equations. Here are the fluctuations in due directly to the noise 

in the equations, and present with zero coupling tOnm = 0, whereas 5^^'^^^^'-^ z^ are the 
remaining fluctuations induced by coupling between modes. The short time assumption 
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was var[log|2;|] <^ 1, which is approximately Kt cosh 2g" ^ 2, from the formal single 
mode solutions (31) A calculation that takes into account nonzero g" is found to give 
the same expression 

The single mode analysis is expected to be accurate only when this ratio 7?. -C 1. 
For this reason, the simulation time improvements due to the local diffusion gauge ()42j] 
are only to be expected while IZ < 1, or perhaps even 7^ <C 1. Let us see what this 
translates to for the case for a uniform gas. 

For a lattice derived from a continuum model, it will be convenient to write the 
lattice interaction strength n = g/hAV in terms of the healing length 

^heal ^ (44) 

This is the minimum length scale over which a local density inhomogeneity in a Bose 
condensate wavefunction can be in balance with the quantum pressure due to kinetic 
effects jSS]- 

In terms of ^^^^^^ the expected simulation time with diffusion gauges is 



^sim ~ O (10)mV^(^^^^^)V^ from §7^, where n is the mean mode occupation. The 
expression can then be evaluated at tgjj^. Imposing 7^ <^ 1 for times t < tgj^^ 
leads to the condition Ax > ^heal^i/4Q ^^-^ -where the quantity Ax = {AVY^'^ is the 
geometric mean of the lattice spacing. Local gauges only lead to significant simulation 
time improvements in the single mode when n ^ 1. 

We conclude that local diffusion gauges will lead to simulation time 
improvements when 

^>0(^heal)_ (45) 

The calculations of Section |S1 will be seen to be consistent with this. This is a strong 
limitation on local diffusion gauges, since interesting dynamical effects can readily occur 
over shorter length scales than this. It suggests that nonlocal diffusion gauges may be 
more useful for general spatially-varying problems, but these are outside the scope of 
the present article. 



4. Local drift gauges 

A complementary approach to the diffusion gauge is to remove the source of the drift 
instability in Eqn. ()225|) by using drift gauges that directly alter the unstable drift 
equations. Since the nonlinearity takes the local form da^ = —inn^a^ dt + ... etc., this 
can be done using drift gauges Qk dependent only on local parameters. These methods 
are certainly capable of eliminating movable singularities, which are a common cause 
of systematic boundary term errors. However, this is accompanied by a corresponding 
growth in sampling error, due to the introduction of a weight term in the stochastic 
equations. 
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4-1- Single mode 

For a single mode with imaginary diffusion gauge as in Section 13 the Ito equations are 

da = a[—iKn — 7/2 — iViK{Qi cosh (7" + iQ2 sinhg") ] dt + i^/ina [dWi cosh (7" + idW2 sinhg"] 

djd = (3[iKn — 7/2 — ViK(—iQi sinh (7" + Q2 cosh (7") ] dt + Vii^P [—idWi sinh (7" + dW2 cosh 
2 

dQ = Qj2^kdWk. (46) 

k=l 

4-1.1. Drift gauge When Qj = 0, the rapid increase in phase variable variance is due to 
the process outlined in point 3 of Section 1231 A drift gauge that interrupts this process 
by removing the offending drift terms dlog |a| = nn" dt and dlog = —nn" dt is 



Qi = iQ2 = —vii^e ^ n" . (47) 

4.1.2. Logarithmic variances There is a price paid, of course, and this is fluctuations 
in the global weight fi, such that var[log|f2|] now scales oc e~^^". With a fluctuating 
weight, quantum phase correlations such as G*-^-* are now given by 

G«(0,t) = /3o(^^«(t))s = al{mt))t (48) 

using (fTUI). Note that since from ()llc|) by the properties of Ito calculus d{Q)s = 0. Thus 
the denominator (i7 + i7*)s appearing in (jlUj) can be ignored if the initial distribution 
is normalized so that (i7)s = 1. From ()48|) . the relevant random variables in phase- 
dependent calculations are (fla) and (fip), so that the log-phase quantity analogous to 
V in the calculations of Section El is now 

Vn = ^(var [log |^^(t)a(t)|] + var [log m)P{t)\]) (49) 

(compare to pO|l ). Proceeding as before by formally solving (j^6|l we find that for off- 
diagonal coherent state initial conditions, 

V. = ^ |t cosh(2/) - e-^>oP ( ^-^1 - e-'^"[{n',r - «)^] ( ^-^^) ] -(50) 



2 [ ' " ' \ q J ^ J ^ 27 

Just as in the pure diffusion gauge case, there is a trade-off between fluctuations 
due to direct phase variable noise oc cosh(2(7") (first term), and fluctuations in the global 
weight n, which are dependent on e~^^" (following terms). The optimum g" = g'^-^^ can 
be calculated by solving dVn/dg" = 0. 

4.1.3. Approximate optimum diffusion gauge Since the aim is a gauge adaptable to 
changing mode occupation like (j^2|) . it is desirable to obtain an approximate expression 
for g'^-p^ that can be rapidly evaluated at each time step in a calculation. To do this, 
we first consider some special cases: 

At short times such that lo'l^opt ^ ^ '^'^^opt^~'^^" ^ -'- optimum is given 

by the roots of the cubic in = e : 

4/ctopt|no|V/ + a3(no,7topt)^/ "1=0 (51) 
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where 

nJnn.v) = 1 +A(n'^)^ ( — 

2v 



asK,!;) = 1 + 4K)2 — -2|r2or • (52) 



For zero damping 03(720, 0) = 1 + 4(nQ; . 

In the case of simulations with sizable occupation of modes, and times up to O {t^^Yi) 
the short time conditions are satisfied and this cubic applies. When the Vg term is 
negligible we have 

9oY>t ~ ^l°g (l^oliy4Ktopt) . (53) 

This occurs when hq is large enough and mostly real, and t is big enough: i.e. when 
'^3('^o,7^) ^ (4/^ioptl^on^''^- For example, when undamped, /t^opt i^^st be at least 
^ 1/4 1 nop with classical initial conditions (rig — 0)- The opposite case when no is 
either too small, too imaginary, or the time is too short has the Vg term negligible and 
leads to 

^opt ~ \ log "3 ("-0 , T^opt ) • (54) 

For strong damping g > 0, we again have linear growth Vn = ntcos\i{2g") /2 — 62, 
where the constant is now 

h, = '-{[{n',r - (n'o')^]e-^^" - \no\V{e'^" - e)]. (55) 

An approximation for the diffusion gauge that is found to work well in practice (see 
Section is 

^approx = I log {4|%P/tiopt + a3('^o, 7^opt)^^^} ' (56) 

which reduces to (I53|) and (jMj) in their limits of applicability. The discrepancy A 
between and the exact optimization obtained by solving dVn/dg" = is shown for 
real no in Figure |21 

4-1 .4- Useful simulation time Consider the undamped high occupation regime with 
coherent state initial conditions no = ^ 1. Using (jS2I), at target times topt ^ 1/4/t^O' 
one has e approx ^ l/[4noKtopt)^^^ ^ 1- Then the terms in ()5U|) of highest order in 
no give Vn(t = ^opt) ~ |[('t''^opt)^^o/4]^^^. (Target times of interest will almost always 
satisfy the prior condition.) We can use the condition Vn < 10 from (PT|) to estimate 
the useful simulation time in this regime as 

W ^ O (lO)tcoh- (57) 
This is again a large improvement compared to the positive P results in Table ^ 



4-2. Extension to many modes 

We now wish to consider how the drift gauge approach can be used to treat a multi-mode 
situation. 
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■■ 0.01 
■■ 0.02 

: 0.05 
: 0.10 



: 0.40 



Figure 2. Discrepancy A = 9 opt ~ 5approx between .g^p^ (the exact optimization 
of g" from (|50|) and the approximate expression (|56|l for an undamped mode with 
diagonal coherent initial occupation no = ti'q. Discrepancy values A are shown as 
contours. Dotted line approximates region of greatest discrepancy Kt^-pt ~ 1/4|'t-o 



4-2.1. Adaptive gauge Proceeding to as in Section ESI the suggested gauges with the 
present approach are: 

g'; = ^log{4|n„(t)p/€trem + a3K(t),7trem)'/'} , (58) 
appearing in noise matrices of the form (jl^, and 

Gj = iQj+M = -ViKlm[nn^ (t)] exp ( -g'' {t) ) . (59) 



for j 



1,...,M. 



4.2.2. Drift gauges and weight spread Drift gauged simulations using encounter a 
scahng problem in many-mode systems because the single weight variable fl accumulates 
fluctuations from all modes (see (jllc|) ). There are precisely two independent noises and 
two drift gauges for each mode. Consider for example, a uniform gas of density p, 
and volume on M modes, so that the mean mode occupation is rT = pAV Writing 
Gk = Q'k + iG'k ^"^^ logfi = 6' + i9" , from (jllcj) one can show that 



d 



-var [6'] = + covar[0', [Qlf] - covar[0', {g',)']),. 



(60) 



For the uniform gas the contribution from each mode is identical on average, and from 
dSD), Q'k and Q'l are of the form ±a/ k/2 e~^"lm [n]. Hence, approximately. 



^var [log \nW oc MK{e~^3"\m [nf)^ . 



(61) 

From the formal solution ()Hla|) (which also applies with drift gauges ()59j) ). one 

has var [Im [n]] = ((|n| sinim [logn])^)s ~ n^(sin^ Im [logn(t)])s, which (at short times 



IS 



d 
~dt 



Kte . Thus at these short times 



var [log f« oc Mn^nhe-^^" 



(62) 
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and at high lattice occupations n ^ 1 when using (|58p e ~ l/(4n/tt)-^/^, so, 
substituting: 

var [log \n\] ^ oc M{gptf'^. (63) 

Imposing the log variance limit 1)211) on \VL\ (because the factor appears in all 
observable estimates ()10|) ). one obtains the estimate that 

1 

^sim gpMV^- ^^^^ 

For this reason we expect that simulations using the local diffusion and local drift 
gauges ()3^ and ()3H|l will only give significant simulation time improvements when the 
number of highly occupied modes is relatively small. Indeed, the two mode cases of 
Section [7| show strong simulation time improvement with this gauge, while we have 
found that the many-mode uniform gases of Section |H1 do not simulate well with the 
gauge form (jSHl)- 

Note, however, that since single-mode tgj^^ decreases rapidly with occupation n, 
it is the most highly occupied modes that limit the simulation time. So, even a very 
large M system may still experience improvements in simulation time with the present 
method if there are only a few modes with the highest occupations. 

As in the pure diffusion gauge case, this does not preclude that better scaling 
may be obtainable with appropriately tailored nonlocal drift and diffusion gauges. In 
particular, the local drift gauge employed here does not take into account the fact that 
neighbouring lattice points with spacings of less than a healing length become strongly 
correlated due to the kinetic energy terms. These cause particle exchange, and hence 
effectively average out local fluctuations in n" . 

5. Single mode: numerical results 

Simulations of an undamped single-mode anharmonic oscillator (see Section 12.3)1 were 
performed for a wide range of initial coherent states uq = nQ from 10~^ to 10^°. Results 
for the standard positive P method were reported in ^2]- Here the gauged methods 
described in Sections El and El are tested. 

5.1. Procedure 

In what follows, the term useful precision for an observable O has been taken to indicate 
the situation where the estimate O = (/)s of (O) using S = 10^ trajectories has a relative 
precision of at least 0.1 at the one sigma confidence level. This is assuming the CLT 
holds so that 



AO^f-^ (65) 

is used to assess uncertainty in O. For the model here, we consider useful precision 
in the magnitude of phase- dependent correlations |G*^^^(0, )f:)|, which is the low-order 
observable most sensitive to the numerical instabilities in the equations. 
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Table 1. Some properties of simulation times tgjj^^ for various gauges when applied to 
the undamped one-mode anharmonic oscillator. Initial coherent state mean occupation 
no = ti'q. 



Drift gauge Diffusion gauge Useful simulation time tsim Maximum uq 
Qk g" maximized over <opt, for which 

when no = tt-q > 1 X^sim > 1 








( 1.06 ±0.16 


) ^osc 


0.014 


+ 0.016 
0.008 







( 1.7 ± 0.4 


^osc 


0.08 


+ 0.07 
0.05 








( 2.5 ± 0.2 




0.11 


± 0.06 





or (36) 


( 8.2 ± 0.4 


^coh 


12 


± 3 







( 10.4 ± 0.7 


^coh 


19 


± 4 


(El 


lESJ 


(30 ±3 


^coh 


150 


± 40 


(EH) 


(t58ll 


(35 ±4 


^coh 


240 


± 70 



Uncertainties in the calculated useful times tgj^^ arise because the AIG*^^)! were 
themselves estimated from finite ensembles of iS = 10*^ trajectories. The range of t^^^ 
indicated in Figure El was obtained from 10 independent runs with identical parameters. 

Taking the analytic scalings (jn7|) (jFfj) at high rio, and from Ref. at low no into 
account, parameters in an approximate curve 

^est = ^|hK)--]-'^^ + 

have been fitted to the empirical data, as was done for the standard positive P method in 
[12]. Best values of Cj are given in Table|21 The expression (jM)) reduces to Ciu'q'''^ / n and 
(cs — C4 lognQ)/(/t/2), when n'^^l and n'^ <^ 1, respectively. Cq determines the high rig 
scaling (here, assumed from analysis of V, or Vq), ci characterizes the pre-factor for high 
n'g, C3 is proportional to a constant residual tgj^^ at near vacuum, C4 characterizes the 
curvature at small Wq, while C2 is related to the stiffness of the transition between the two 
regimes. Uncertainty Acj in parameters Cj was worked out by requiring ^^^{[testl'^i =t 
Acj,no) - WK)]/Atsim}^ = E„o{l + ([W(ci,no) - tsim("'o)]/^W)^}- 

For calculations involving a diffusion gauge dependent on target time ^opt' ^ wide 
variety of target times were tried to investigate the dependence between t^^^ and ^opt' 
and ascertain what are the longest simulation times achievable. 



2 log 



/4q 



-C2 



(66) 



5.2. Simulation times 

Figure 01 compares tgjj^, as defined via "useful precision" in |GW(0,t)|, for several gauge 
choices. Note that a logarithmic scale is used. Results at high occupation are tabulated 
in Tabled which includes data for a larger set of gauge choices. Figure |3 gives examples 
of calculated values |G'^^^(0, t) \ along with error estimates. Table|21gives empirical fitting 
parameters to the expression (jHEl). 
We see that: 
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Figure 3. Maximum useful simulation time igjjjj, of the one-mode undamped 
anharmonic oscillator with various gauge choices. Initial coherent state mean mode 
occupations no = "-o- Width of plotted lines shows range, when using 10 runs of 
S = 10** trajectories each. The drift gauge is l|59|l . while the diffusion gauge is H42(l 
when on its own, or H58|l when combined with the drift gauge. 

Table 2. Empirical fitting parameters for maximum useful simulation time igjj-^ with 
several different gauge choices when applied to the undamped one-mode anharmonic 
oscillator. The fit is to expression (|56|l . 





Positive P 


Drift gauge 


Diffusion gauge 


Both gauges 


Gfc 





(EH) 





(El 


g" 








(EH) 




Co 


2/3 


1 


1/2 


1/2 


Cl 


2.5 ± 0.2 


11 ±3 


10.4 ± 0.7 


35 ±4 


C2 


3.2 ±?°2 


T A _|_ oo 
=■= 0.4 


9 7 _(_ oo 
=■= 1.0 


3.6 ±^3 


C3 


-0.5 ± 0.3 


-0.5 ± 0.3 


2.4 ± 0.6 


2.8 ± 0.9 


C4 


0.45 ± 0.07 


0.49 ± 0.08 


0.23 ± 0.13 


0.23 ± 0.13 



• Combining drift and diffusion gauges gives the longest useful simulation times. Such 
simulations give good precision well beyond the point at which all coherence has 
decayed away for highly-occupied modes - potentially up to about 35 collapse times 
tcoh in the cases treated here. 

• Diffusion-gauge-only simulations also give quite good statistical behaviour 
(although useful simulation times are about 4 times shorter at high occupation 
than with both gauges). 

• Despite the efficient behaviour of combined gauge simulations, using only a drift- 
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Figure 4. Modulus of the phase correlation function G'^^^(0, i). Comparison of 
calculations with different gauges. Subplots (b) and (e): use diffusion gauge H68|l of 
Plimakei al [21] with t^-p^ = '^^coh^ while subplots (c) and (f)use the combined drift 
and adaptive diffusion gauges (|^ and with the choice tgpt = '^^^coh- "^^^ initial 
conditions: coherent state with (n) = rip. Triple solid lines indicate G^^^ estimate with 
error bars, dashed lines are exact values. S = 10"* trajectories in all cases. 



gauge gives even worse statistical error tlian no gauge at all. Such simulations are 
restricted in time to about one phase oscillation. This indicates that the drift gauge 
choice made here is not an optimum choice when used in isolation. 

• At low occupation, i.e. of the order of one boson or less, combined gauge methods 
still give the best results, but the advantage is marginal. 

• The simulation times with nonzero diffusion gauges (whether accompanied by drift 
gauge (j59|l or not) not only have better scaling with no when uq is large, but this 
power-law decay of simulation time starts much later, as seen in Figure El and the 
right hand column of Table ^ 

5.3. Target time dependence 

Figures El and El show the dependence of simulation times on the target time parameter 
^opt ^ variety of gauges. Comparison is made between the full adaptive forms 
g"{nn(t),trem) (solid lines), number-adaptive-only forms 5'"(^n(i^), i^opt) ^^^^ 
use explicit time-dependence trem (dashed lines), and the Plimakei al jl^ gauge 
(dash-dotted lines). 
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^opt/ '^coh topt/tcoh 
Figure 5. Dependence of actual useful simulation time t^^^ on the a priori target 
time tgpt ^'^^ ^ variety of diffusion gauges in the Gk = schemes: Solid line: full 
adaptive gauge l|42(l : dashed line: number-adaptive-only gauge of form (36) but with 
the modification no n{t); dash-dotted line: gauge (|68|l of Plimakei al For the 
gauge (36), the whole region where useful precision occurs is shown by the dotted line. 
Undamped anharmonic oscillator system (|13|l . with coherent state initial conditions 
no = nL 
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Figure 6. Dependence of t^-^^ on tppt' ^ '^"^ Figure |31 with drift gauge as well 
as diffusion gauges. Details as in Figure El 



Some comments: 

• The diffusion gauge forms that optimize for the "remaining time to target" trem 
give the longest simulation times, and these times are well controlled. Statistical 
error can be reliably expected to remain small up to the exphcit target time ^opt' 
provided that this is within the useful simulation range given in Table ^ 

• Diffusion gauges that instead use a constant tQp^ parameter give: 

(i) Somewhat shorter simulation times. 
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(ii) A complicated relationship between target and useful simulation times. 
Broadly speaking ^Qp^ <C 'f^sim optimum cases. 

These forms of diffusion gauges require tedious parameter searching to find the best 
t^-pf^ choice for given initial conditions. The probable reason why t^^^ ^ ^Qp^ here 
is that g" has been optimized to minimize the variance of logarithmic variables. 
This is not the same as the variance of the non- logarithmic variables a, or f3 that 
actually appear in the observable calculations. Hence, different g" forms may extend 
simulation time by a further amount. 

• When drift gauges are used, an adaptive diffusion gauge g"{n{t)) rather than a 
constant g"{nQ) can give much longer simulation times. (Compare to the Plimak et 
al gauge in Figure ini(a) which effectively scans through a range of constant g" 
values) . 

• At low nQ, tgjj^ is only weakly dependent on local diffusion gauge choice. 

• When there is no drift gauge, for tiq ^ 1, the time-adaptive gauge forms g"{tiem) 
lead to a peculiar effect if the optimization time tgp^ is chosen larger than the usual 
maximum tgjjj^ given in Tabled The statistical error in the G*-^^ estimate first rises 
rapidly, then falls again, and finally grows definitively. The parameter region in 
which this occurs is shown in Figure EK a). In effect one has two time intervals 
when the simulation gives useful results: at short times, and later in a time interval 
around t ~ O (lOtp^j^). 

More detail on these numerical investigations can be found in Ref. ^Hj, Chapter 7. 
5.4- Comparison to recent related work 

Improvements to the basic positive P simulation method for specific cases of interacting 
Bose gas systems have been tried with some success in several recent publications 0123 
OH Ej . Here we compare these with the stochastic gauge formalism, and make some 
comparison to the results and analysis in the present chapter. 

5.^.i. The work |^ of Carusotto, Castin, and Dalibard An isolated (i.e. particle- 
conserving) system of exactly N interacting bosons was considered (on a ID lattice). 
The "coherent state simple scheme" for stochastic wavefunctions described in Section 
III B 2 therein can be identified as using drift gauges of the form 

when re- written our notation for the single-mode system (jl3|) . This gauge causes a 
full decoupling of the complementary a and /3 equations by making the replacements 
n— i> |ap orn— i> in the nonlinear terms. Like ()59|) it is also successful in removing 
movable singularities, since the nonlinear terms in the radial equations for and d\j3\ 
are removed. 
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The stochastic wavefunction in the form 
systems with a definite and exphcit number 
coherent out-couphng from a system cannot 



presented in j8j is only apphcable to closed 
of particles, so e.g. evaporative cooling or 
be treated. 



5.4-2. The work of Plimak, Olsen, and Collett A single-mode undamped, 

gainless system (fT^ at high Bose occupation with coherent state initial conditions 
was considered. The "noise optimization" scheme applied therein to greatly improve 
simulation times can be identified as an imaginary diffusion gauge of the form (rewritten 
in the present notation) 

g" = g\ = ]^ cosh"^ ^o«;i^opt (68) 

defined at high occupation or long times (i.e. while ?^o'^^opt — This is dependent on 
a target time ^Qp^ (which was taken to be ^Qp^ = 3^^^]^ in the calculations of Ref. j^). 
and the initial Bose occupation no = u'q. The useful simulation times obtainable with 
this method are also shown in Figures El and El and Table (H Their dependence on t^-^i 
has been calculated here, and is shown in Figure El 



5.4.3. The work /?/ of Drummond and Deuar In section 5.3 of the above article, some 
preliminary results for the dynamics of a one-mode, undamped, gainless system (with 
no = 9 particles on average) were shown. The drift gauge ()59|). and a constant imaginary 
diffusion gauge g" = 1.4 were used. 



6. Convergence issues 
6.1. General 

A subtle issue with many phase-space distributions is the possibility of so-called 
boundary term errors. These can arise when the tails of a distribution (say G{a, (3, Q)) 
do not fall off fast enough as the boundaries of phase space are approached (in the 
case here, as or \Q\ — > 00). It is possible for this to lead to a bias in means 

of random variables even in the infinite sample limit, if parts of the distribution which 
have a non-negligible effect are never sampled. 

Some numerical indicators have been developed j3ti] that allow one to search for 
symptoms of these errors using the numerical data. The most useful of these indicators 
is sudden appearance of spiking in observable estimates, where the onset of this spiking 
tends to come earlier as more trajectories are added. Such spikes usually occur when a 
single trajectory samples a long power-law tail. Results obtained after first spiking 
are suspect. 

In our experience [1,6|, another indicator can be obtained by performing two 
simulations with sample sizes S differing by an order of magnitude. If a statistically 
significant difference (e.g. 2cr) in observable predictions occurs, the simulation 
is suspect. 
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In all cases with which we are familiar, symptoms of potential convergence problems 
have been apparent already in the equations of motion. They have been of two 
kinds TFjf: 

6.2. Divergence symptoms of the first kind: movable singularities 

This is a symptom apparent in the drift parts of the equations of motion. 

It has been found jHE] that a systematic boundary term error is often associated 
with the presence of so-called movable singularities. These are trajectories (usually of 
measure zero) in the deterministic parts of the equations that diverge in a finite time. 
The presence of a bias in this situation can be understood by considering that the effect 
of an infinitesimal number of divergent trajectories may be nonzero. Some cases where 
this occurs with a positive P simulation have been investigated by Gilchrist et al |36j . 
The bias has been shown to disappear when the movable singularities are removed using 
drift gauges jni- See also Ref. jTHj, Chapter 6. 

Consider a set of generic stochastic equations 

dv = A,{^r)dt + J2 B,k{^r)dWk{t), (69) 

k 

in complex variables v = {v} free to explore the whole complex plane. If, in the limit 
|f I — ^ oo, deterministic growth of \v\ is exponential or slower, then the trajectory cannot 
reach infinity in a finite time by deterministic processes, and moving singularities are 
ruled out. Exponential growth in v occurs when Re [A^/v] is a positive constant, so we 
conclude that a condition sufficient to rule out moving singularities is that 

lim Re — (70) 

converges for all variables v. 

6.3. Divergence symptoms of the second kind: noise-weight divergence 

This is a symptom which arises from a combination of the noise behaviour and the form 
of the quantities which are averaged to obtain observable estimates. 

A classic and well known example of this convergence problem occurs for the present 
single mode system ()13p when using an un-normalized Bargmann coherent state kernel 

A(a,/5) = (71) 

(Compare to the gauge P kernel (jSl) )• The (non) convergence of this model has been 
considered by Carusotto and Castin|2ZI, and also in Ref. [TE], Section 6.2.2. 
One finds that the Ito equations of motion are 

da = ia\/iK dWi (72a) 
dp = pVii^dW2. {72b) 

f Of course, other ways of achieving a divergence may occur 
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To calculate estimates for an observable O we need to evaluate 

(73) 



Tr 


'dp 




AO 




Tr[p] 




A 





for some appropriate /(a, /3) which depends on the details of O. The equations of motion 
can be formally solved to give 

n{t) = a{t)l3{t) = n(0)exp [-^S,- {t) + i^e{t)\ , (74) 
with the variance t Gaussian random variables defined similarly to ()31c|) : 



e = ^f [dW,{s) ± dW,{s)] . (75) 
V 2, Js=o 

For any observable estimate, we must estimate (e"^)s using our samples of n{^^,^^). 
The distribution of this can be explicitly evaluated: 

/oo 
P(e+)P(r) e"f'd^^d^2 
-oo 

oc 1^ exp |n(0) e-^'^'e^^'^^ " ^ " ^} ^^^^ 

This distribution is divergent as — — oo because the e^^^ factor in the 
exponent always beats the (relatively) weak gaussian tail, which can only manage a 
quadratic — (^~)^/2t drop-off in the exponent. 

6.4- Stochastic model system 

A more general feeling for the effect on convergence of the noise-weight relationship can 
be gained by considering the following simple equation for a complex variable v. 

dv{t) =cv{trdW{t), (77) 

with real constants n and c. Its formal solution is 

/ b(0)^"" + c(l-n)e(t)]^ ifn^l, .^g. 
\ v{0) exp[c^{t)] ifn = l. ^ ^ 

where again ^(t) = Jg^QdW{s) is a Gaussian random variable of variance t, mean zero. 
One finds that: 

, ^, [ converges Vm and Vt for n < 1 
I does not converge for > 1 

and 

converges Vm < 2(1 — n) and Vt for n < 1 

(exp(f™'))s ^ converges for m = 2(1 — n) and for t < 2c{i-n)^ for n < 1 (80) 

does not converge for n > 1 
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One can see from the above that, barring special favourable circumstancesf, the 
following relationships will hold (the notation is as in (jUIIj) ): 

A) If observable averages involve only expressions polynomial in the variables v, 
then schemes for which noise terms grow faster than linearly as |f | — oo will be 
divergent. That is those where 

lim Re[S„(v)/t;] (81) 

is unbounded for any v. For example = cf where n > 1. 

B) If observable averages involve expressions exponential in variables v, then 
schemes for which noise terms grow faster than oc \/\v\ as — > oo will be divergent. 
That is those where 

lim |5^(v)/v^| ^ 0. (82) 

The special case when the limit (I82p is finite is also divergent, but only after some initial 
time period. 

6.5. Single mode: diffusion gauge 

How do the gauged methods compare to these divergence symptoms for the single-mode 
system? The diffusion-gauged method of Section El has a finite formal solution (31) at 
all times. For observables of finite order in and a, only polynomial expressions of a 
and P need be averaged. Such expressions scale as oc exp [factors x C,j(t)], and thus their 
stochastic averages are equivalent to integrals like J P{C) exp[factors x ^] dC, which are 
convergent due to the Gaussian form of -P(0- 

6.6. Single mode: drift gauge 

The equations ()46|) with the drift gauge ()47|) can be formally solved: 

logn(t) = logn(O) - 7t - e'^" ^/^ [Cit) - z^+(t)] (83) 

logn(t) = - e-^"v^ [ Im [n{t)] [dWi{s) - tdW2{s)] . (84) 

Js=0 

(taking Q{0) = 1). From (fTUI) . observables require the averaging of quantities like 
Qf{a,f3), which include the factor |^^(t)|. This takes the form 

\Q{t)\ = exp j-v^e-^'" ^ |n(0)|e-^^e-^""'"«"(^) sin [Zn(0) + v^t{s)] c/W^+(s)| , (85) 

where dW+ = {dWi + dW2)/V2. 

For large negative ^~(s), the gaussian drop-off of P(^^(s)) is insufficient by itself 
to directly prevent the divergence of this factor. However, the situation is more subtle 
than in the simple model equations, due to the presence of an oscillatory stochastic 

f e.g. a topological barrier which prevents any trajectories reaching \v\ — s- oo, given the right initial 
conditions. 
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term in the integral. Hence the convergence of the observable averages is still an open 
question for the drift-gauged method, despite the absence of movable singularities. 

Lack of convergence has been shown for the "simple coherent" stochastic 
wavefunction method|SZl, which has been discussed in Section 15.4.11 



6.7. Many-modes 

Regarding movable singularities in the many- mode equations (11), let us first compare 
to the movable singularity condition (fTOj) . One sees that: 

1) The drift terms dependent on Unn and 7n as well as the noise terms dependent 
on K lead to only exponential growth in the moduli of and and so can not cause 
movable singularities. 

2) The mode-mixing terms dependent on Unm, where n 7^ m, lead to drift of the 

form 

d\an\ = |am||t^nm| siu (ZcUnm + ^ttm - ^ttn) dt + . . . 

< |«m||t^nmMt + . . . . (86) 

and of similar form for |/5n|- These terms lead to behaviour like a linear matrix 
differential equation for the radial evolution of the coherent state amplitudes. Such 
equations do not diverge in finite time, their solutions being finite linear combinations 
of exponentials, hence these terms by themselves can not lead to movable singularities 
either. 

3) Only the nonlinear drift terms can cause movable singularities. With no drift 
gauge {Qk = 0) as in the positive P or diffusion-gauge-only equations, these lead (when 
dominant) to evolution of the form 

d\an\ = K\an\^\Pn\ sin(Zan + Z/3n) + • • • 

d\Pn\ = -i^lPnl |an| sin(Za;n + Z/3n) + . . . . 

These can cause either or |/3n| to diverge at finite time, for some trajectories. For 
example, if |/5n| sin(ZQ;n + Z/?n) = K conspire to be approximately constant and nonzero 
over the relevant timescale, then 

1 - KK\an{to)\{t - to) 

which diverges at time tdiv = + l/'^-^l'^n(^o)|- Of course the precise condition ''K 
is constant over the time to to t^iv" will only occur for a set of trajectories of measure 
zero, but this is typical of movable singularities. 

4) The equations that use the drift gauges (|5^ do not contain the terms (|H7j). and 
no movable singularities occur. 

As for noise-weight divergences, the situation appears to be similar to the single- 
mode cases. The formal solution for the weight is 

f2(t) = exp < —iViK J ^^Im [^^n(■s)] x noise increments(s) > , (89) 
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and contains an exponent containing exponentials of gaussian random variables (s) 
as part of the n^is). This is indicative of possible divergences. 

6.8. The "what happened to the divergence?" puzzle 

It appears that the observable averages using the drift gauge (jlTjl have the potential to be 
non-convergent even in the single mode case. Why then are the numerical simulations 
of Section well behaved for such a long time, showing no sign of bias? Similarly, 
no systematic error was seen in simulations using the "simple coherent" stochastic 
wavefunction scheme |H], despite the subsequent proof of its divergence in Ref. ^Tj . 

The detailed investigation of this is beyond the scope of this paper. The 
simplest explanation is simply that the distribution tails are sufficiently convergent to 
eliminate boundary terms, while still having a large (perhaps infinite) variance in some 
observables. It is possible that the appearance of large statistical uncertainty masks 
any systematic errors that may occur in the S ^ oo limit. This is plausible because 
long distribution tails certainly give rise to large phase-space excursions and thus huge 
statistical uncertainty, whether or not systematic biases in the limit S ^ oo are present. 
Another reason for this lack of bias may be some type of special symmetry properties 
in the equations or the observables calculated. 

6.9. Summary 

Divergences of the moving singularity type may be present in many-mode (but not single- 
mode) simulations using the diffusion-gauge-only method of Section IHl while divergences 
of the noise-weight type may be present in simulations with the drift-gauged method 
described in Section H] However, no bias of any kind was seen in the single-mode 
simulations (or the two- mode simulations, as shall be seen below). 

Hence, the simulations appear to give correct results, but numerical indicators 
such as spiking and ensemble-size dependence should be rigorously monitored in all 
calculations. 

7. Two coupled modes 

We now look at the behaviour of a two-mode system, to investigate how the statistical 
behaviour seen for the single mode is affected once coupling between modes is present. 
This is a simple enough system that investigation of several examples gives meaningful 
insight into the general situation. 

7.1. The model 

The system consists of two orthogonal modes labeled 1 and 2 with inter-particle 
interactions in each mode, and Rabi coupling between them. No damping will be 
considered, for simplicity. The coupling frequency will be restricted to be real. 
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Time units are chosen so that the nonhnear interaction frequency is k = 2, and a 
transformation is made to an interaction picture in which the hnear mode self-energies 
hjjjjjUj are moved into the Heisenberg evolution of operators. The interaction picture 
Hamiltonian then is 

2 

H = fVjJi2 

i=i 

The Rabi frequency is uji2 (in scaled time units). Stochastic equations are as in ()2.2|) . 
Physically this model can represent, for example, two internal boson states coupled 
by an EM field, or two trapped condensates spatially separated by a barrier. This 
approximation has been widely used to investigate the quantum behavior of BECs in 
two-state svstemspS]. 

Let us consider two kinds of initial conditions which broadly represent the two kinds 
of situations generically occurring in all many-mode simulations: 

(i) Case 1: Coupling between modes of widely differing occupation. 

(ii) Case 2: Coupling between modes of similar occupation. 

In a many- mode calculation, adjacent modes typically behave like case 2, since if a field 
model is well resolved by the lattice, then physical properties (e.g. density, and hence 
mode occupation) should not change much over the distance between neighboring lattice 
points. Long distance coupling will tend to behave like case 1. 



(90) 



7.2. Case 1: Coupling to a vacuum mode 

The system starts initially with a coherent state of mean particle number in mode 1, 
and vacuum in mode 2. In all simulations of this case, the inter-mode coupling strength 
was taken to be uji2 = 5, but the mean particle number uq (conserved in time) was 
varied. 

At low particle number, the Rabi oscillations dominate the Hamiltonian, and 
particles oscillate between the modes, without much phase collapse. At high particle 
number uq, on the other hand, phase collapse dominates mode 1, suppressing also the 
coherent transfer of particles to mode 2. 

The particular values chosen to simulate were 

no = {1,17,200,1500,10^}, (91) 

Simulation times were assessed using the calculated uncertainties in the two observables: 
1) The fraction of particles in the (initially empty) mode 2: 

P. = (92) 
where rij = OjUj, and 2) the local normalized second order correlation functions: 
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The second order correlations quantify the amount of (instantaneous) bunch- 
ing/antibunching in the boson field. These are unity for coherent states, two for ther- 
mal fields, and 1 — 1/n for Fock number states of n particles. Large values g^"^^ > 2 
occur e.g. for quantum superpositions of vacuum and Fock number states with two or 
more particles where the average particle number is small. For example in the state 
= sin^lO) +cos^|n), ^(2) = (1 _ i)/cos2^. 

The drift and diffusion gauge scheme using and ()58|) was considered. 

Comparison was also made to the positive P {g" = Qk=0), and to the special case 
of tQp^ = in the diffusion gauge which then is 

9: = l^og[l + Anl{tr]. (94) 

This may become nonzero after spread in the from the initial n" = occurs. In each 
run, iS = 2 X 10^ trajectories were used, and useful simulation precision taken to occur 
at such a time tgim when 10% or smaller relative uncertainty in an observable could be 
obtained using 10® trajectories. 

In Figure [3 simulation times are compared to physical timescales and expected 
values based on single-mode expressions from Table El An example simulation is 
shown in Figure |H1 Note that for gauged simulations using and (jHH), single-mode 
simulations led to the tgjj^ empirical fitting parameters 

{co, . . . , C4} = {1, 800 ± 260, 3.6 ^^3, 1.2 ± 0.2, 0.42 ± 0.03}. (95) 

One point to note is that t^^^ was based on the moment when relative error in a 
quantity gj or pj was first found to be too large. In calculations of g^'^\ uncertainties 
are much greater when g^'^\0, t) peaks — see e.g. Figure |HKf). Good accuracy can often 
be obtained between peaks for much longer times than shown in Figure [TKc) or (d), up 
to about the same simulation time as worked out based on p2. This is especially evident 
for riQ = 10"^. 

7.3. Case 2: coherent mixing of two identical modes 

The system starts initially with identical coherent states of mean particle number hq in 
both modes. The two-mode state is separable. This time simulations were carried out 
with constant particle number uq = 100, but the coupling frequency uj^ was varied. 

At low frequency <^ = 100, phase collapse local in each mode dominates, and 
phase oscillations in each mode occur with period tosc = Tr/riO) while at high frequency 
UJ12 nQ = 100, the inter-mode coupling dominates and phase oscillations for each 
mode occur with period t-^^^^ = It^Iojyi- One expects that for weak coupling the two 
modes should behave largely as two independent single modes of Sectional 

The particular values chosen to simulate were 

u;i2 = {5000, 500, 50, 5, 0.5, 0.05, 0.005, 0.0005}. (96) 

The simulation time t^^^ was assessed here based on 10% uncertainty in |G''^^)(0, t)| 
(for either mode), as in the single- mode case. There were S = 10^ trajectories per 
simulation. 
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Useful simulation times t^YTn f^'" ^ mode coupled to vacuum as in 
Section FT^ Calculated simulation times are shown as data points, with the symbols 
denoting gauge used. positive P; "0"^ drift gauge fKi^l and diffusion gauge 

()94|l : "v"- drift and diffusion gauges (|59|l and H58|l with best target time parameters: 
^i2^opt = {2.5,0.5,0.5,0.075} for no = {1,17,200,10'*}, respectively. Dependence of 
^sim '^"^ ^opt '^^^ found in 16 , Figure 8.6. Subplots (a) and (b) show simulation 

(2) 

times based on estimates of the observable p2, while (c) and (d) times based on g\ . 
Subplots (a) and (c) compare to physical time scales, including Rabi oscillation period 
^Rabi ~ 27r/a;i2, while subplots (b) and (d) compare to expected simulation times 



for a single mode using the empirical fits of Table |2 and H95|l . The expected t^Ym. 
are plotted as light lines. Dotted: positive P; solid: drift and diffusion gauges with 
'•opt ~ 0; dashed: with optimum topt choice. 



In Figure IHl simulation times are compared to physical timescales and expected 
values based on single- mode expressions using Table |21 An example simulation is shown 
in Figure ITUl 

Note that the t^^^ = gauge performed better than any t^-p^ > gauge for 
all uji2 values tried here. The dependence of tgjj^ on ^Qp^ is shown in Figure 8.8 of 
Ref. P|. 
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LJl2t UJut iOut 

Figure 8. Coupling to vacuum mode as in Section FT^ no = 10"'. Subfigures (a)-(c) 
show results with positive P simulations, whereas (d)- (f ) show results with combined 
drift gauge and diffusion gauge l(5H|l using target time (^i2tQp^ = 0.1. Triple lines 
indicate mean and error bars. 

7.4- Analysis 

The above examples have not by any means been a comprehensive assessment of gauge 
performance for general cases of the model since only a few parameter regimes 

have been explored. Still, several aspects of the situation when modes are coupled have 
been seen: 

• Broadly speaking, when local scattering within a mode dominates over the coupling 
between modes, the response of the system to local gauges is similar to what was 
seen for single modes. Large extensions of tgjj^ are found relative to the positive 
P method to well beyond t^^]^ when n ^ 1. Little improvement is found when 
n < 0(1). In this strong scattering regime, one doesn't need much inter-mode 
coupling uji2 to reduce the simulation time in absolute terms by a factor O (2 — 10), 
although tgjj^ > t^^Yi is still obtained. 
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Figure 9. Useful simulation times tgjj-^ for a two identical modes undergoing coherent 
mixing as in Section 17.31 Calculated simulation times are shown as data points, with 
the symbols denoting gauge used. positive P; "0"^ drift gauge and diffusion 

gauge (|94|l . Subplot (b) compares to expected simulation times for a single mode using 
the empirical fits of Table [3 and H95|l . These expected tgjj-^ are plotted as light lines: 
dashed: positive P; solid: drift gauged. 
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Figure 10. Mixing of two identical modes as in Section l7T51 U!i2 — 0.0005. Subplot (a) 
shows results with a positive P simulation, whereas (b) shows results with combined 
drift gauge (|59|l and diffusion gauge H94() . Triple lines indicate mean and error bars. 



When the inter-mode coupling dominates, the local diffusion gauges do not appear 
to be useful. They actually reduce simulation time as compared to the positive P 
method, although several Rabi oscillation periods can always be simulated. In the 
Case 2 simulations, the transition between the strong and weak coupling behavior 
appears to be at around 

iUu^O{Kn), (97) 

which is the point at which the expectation values of the coupling and two-body 
scattering energies in the Hamiltonian are approximately equal. This implies, as 
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noted previously, that nonlocal gauges are likely to be more useful in these cases. 

• The beneficial effect of choosing tgp^ > seen for a single mode (see e.g. Figure EI) 
appears to be suppressed at intermediate mode occupations ^ O (200). However, 
the t^p^ = diffusion gauge has a marked beneficial effect even when 

two body scattering is significant. This parameter range is for uJi2 <^ nn, and 
O (1) < no < O (10^). Simulation times obtained are smaller by about a factor 
of O (2) than those given for a single mode with the same gauge. At higher 
occupations, the benefit gained with (|MI) abates but nonzero ^Qp^ values appear 
to become useful again, and continue to provide strong improvements over positive 
P simulations. (See, e.g. Figure [TKd).) The gauge (|M|) is convenient also because 
there is no a priori parameter topt- 

At the level of the stochastic equations, the evolution of daj can gather randomness 
from three sources 

(i) Directly from the local noise term oc aj^fKdWj 

(ii) Indirectly from the local nonlinear term oc najUj that can amplify variation in the 
local noise term. 

(iii) From the other mode through the coupling term oc Uud^j- 

The drift gauges (|59|l neutralize source 2. The diffusion gauges (H2|l . (j58|) or (jHlI) 
suppresses the direct noise source 1. No local gauge is good at suppressing the third 
source of randomness, however, because these fluctuations are largely independent of 
any processes occurring in mode j. What happens is that even small randomness in one 
mode feeds into the other, can become amplified, and fed back again. Combating such 
effects would require a nonlocal gauge. 

Lastly, significantly more detail on two-mode simulations, including consideration 
of some other gauges, and the relationship between ^Qp^ and tgjj^ can be found in 
Chapter 8 of Ref. 

8. Many-mode example: uniform gas 

We revisit the uniform gas system simulated in the companion paper This consists 
of a uniform one-dimensional gas of bosons with density p and inter-particle s-wave 
scattering length a^. The lattice is chosen with a spacing Ax ^ as so that inter-particle 
interactions are effectively local at each lattice point, and the lattice Hamiltonian (Q) 
applies. Periodic boundary conditions are assumed. 

The initial state is taken to be a coherent wavefunction, which is a stationary state 
of the ideal gas with no inter-particle interactions (i.e. Os = 0). Subsequent evolution 
is with constant > 0, so that there is a disturbance at t = when inter-particle 
interactions are rapidly turned on. Physically, this can correspond to the disturbance 
created in a BEC by rapidly increasing the scattering length at t ~ by e.g. tuning 
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Figure 11. Simulation times as a function of the target time tQ-p^ when using the 
diffusion gauge (|42|l . and = 0. Standard positive P method when tgp^ ~ 0. All 

simulations are of a p = 100/^^^^^ gas, but with differing lattice spacing Ax. S = 10^ 
trajectories, M = 50. 



the external magnetic field near a Feshbach resonance W)\ . More discussion of the 
physics can be found in the companion paper [T^. 

A brief set of exploratory simulations were made with the diffusion gauge ()42j] 
(no drift gauge) to investigate whether simulation time can be extended. The case 
of p = lOO/^'^^^^ was picked, and target time tgp^ was varied for the two cases 

Ax = ^^^^^/2 and Ax = 10^^^^^. Note that in both cases the occupation per mode 
n = pAx is ^ 1. Simulations were with M = 250 and M = 50 lattice points (meaning 
12500 and 5 x 10"^ particles on average), respectively, and S = 10^ trajectories in both 
cases. 

Simulation times obtained are shown in Figure ^2 and correspond to the time of 
first spiking seen in estimates of the second-order spatial correlation function 

—(2)/ \ _ J_ \ ^ (o-m'^m+nQ-mQ-m+n) /qoN 

The time scale used is 

m(e^eal)2 ^ 

(the "healing time"), which is approximately the time needed for the short-distance 
O (i^^^^^) inter-atomic correlations to equilibrate after the disturbance. See e.g. 
Ref. [ig. 

For comparison, gauged simulations of isolated lattice points with the same 
occupations n = pAx gave much bigger improvements in tgjj^. Using (jUBj) and diffusion 
gauge (021) data from Table El we expect tgj^ ^ 140^^ and t^^^ ^ 660tg for Ax = 2^^^*^^ 
and lO^'^^^^, respectively. This is to be contrasted with the predictions of tgj^^ ~ 15t^ 
and tgjj^ ~ 49tj, respectively, for the standard positive P method, which are also an 
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accurate reflection of tgjjj^ in tlie present many-mode case. Tliis is consistent witli wliat 
was seen in Ref. |12j . 

Clearly, the local gauges extend simulation time in these large uniform density 
systems, although the improvement is relatively smaller than for the single mode. Also, 
the numerical results are consistent with the analysis in Section which indicated in 

that simulation time time improvements would occur once Ax > O (e^eal) . Lastly, 
it is noteworthy that the dependence of tgj^^ on t^-p^ shows qualitative similarities to 
the single-mode case of Figure El 

9. Conclusions 

The positive P representation method is capable of simulating many-body quantum 
dynamics of interacting Bose gases from first principles |14 [ I15 | IT2 | I16j . A limiting factor 
is that precision is lost after a certain time tgj^^ due to sampling error, and this may 
be accompanied by systematic boundary term errors. This time can be short when 
some modes of a many-mode boson system are occupied by many particles. In that 
case, precision is lost before these highly occupied modes lose coherence due to phase 
diffusion [HI. 

Using the related gauge P representation, local diffusion and drift gauges have been 
developed that greatly improve the useful simulation times at high occupation, most 
strongly for single-mode cases. The resulting simulation times have been investigated 
in some detail, and substantial improvement has been demonstrated in 1, 2, and many- 
mode (M = 50 and M = 250) simulations. For many of the single- and double-mode 
cases considered, full decoherence can be achieved while still retaining good precision. 

Two gauge choices were proposed as being the most advantageous: 

(i) Diffusion gauge only. 

(ii) Drift and diffusion gauge (j3H|) . 

The latter appears less broadly applicable because weight fluctuations accumulate from 
all modes, leading to a decrease in simulation time with increased system size, and 
because of some doubt over its convergence properties. However, this latter method can 
lead to longer simulation times when there is one, or only several, dominant modes. 

Considerations of convergence in Section IHl have highlighted two routes that may 
lead to divergences: moving singularities and noise-weight considerations. Nevertheless, 
no sign of any bias is seen, an interesting situation which warrants further investigation. 

Generally speaking, the local-gauge methods considered here give much better 
performance when the inter-particle scattering local to each mode is a stronger process 
than inter- mode coupling due to kinetic evolution and/or external potentials. This is 
the case if sufficiently coarse lattices are used. 

Finally, the investigation of local gauges carried out here may also be useful in 
developing more robust nonlocal gauges. These might lengthen simulation times for 
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lattices denser than the heahng length, beyond what is possible with either the standard 
positive P method, or the local stochastic gauges. 
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